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Abstract. In this paper, we consider the flow of an incompressible fluid in a deformable porous 
sohd. We present a mathematical model using the fi-amework offered by the theory of interacting 
continua. In its most general form, this framework provides a mechanism for capturing multiphase 
flow, deformation, chemical reactions and thermal processes, as well as interactions between the 
various physics in a conveniently implemented fashion. To simplify the presentation of the frame- 
work, results are presented for a particular model than can be seen as an extension of Darcy's 
equation (which assumes that the porous solid is rigid) that takes into account elastic deformation 
of the porous solid. The model also considers the effect of deformation on porosity. We show that 
using this model one can recover identical results as in the framework proposed by Biot and Terza- 
ghi. Some salient features of the framework are as follows: (a) It is a consistent mixture theory 
model, and adheres to the laws and principles of continuum thermodynamics, (b) the model is 
capable of simulating various important phenomena like consolidation and surface subsidence, and 
(c) the model is amenable to several extensions. We also present numerical coupling algorithms to 
obtain coupled flow-deformation response. Several representative numerical examples are presented 
to illustrate the capability of the mathematical model and the performance of the computational 
framework. 



1. MOTIVATION AND INTRODUCTION 

Higher oil prices and increasing awareness of the environmental impact of carbon pollution 
have motivated substantial interest in carbon-dioxide capture and storage (CCS). There is a growing 
consensus in both policy circles and in the energy industry that within the next few years, the US 
federal government will adopt some form of regulation for CO2 emissions [1]. At the same time, it 
is widely believed that much of the energy supply over the coming decades will continue to come 
from fossil fuels. Many analysts believe that the only way to reconcile the anticipated growth in 
the use of fossil fuels with anticipated limits on CO2 emissions is through the development and 
deployment of carbon-dioxide capture and sequestration. 
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Figure 1. Lake Nyos is a deep crater 
lake in Cameroon (top figure). In 
1986, huge amounts of carbon-dioxide 
gas were released from the lake possi- 
bly due to seismic tremors. A mix of 
carbon-dioxide and water erupted 120 
meters above the lake with a speed of 
roughly 100 kilometers per hour. This 
displaced oxygen resulting in asphyxi- 
ation of more than 1,700 people and 
3,000 cattle as far as 23 kilometers 
from the lake (bottom figure). [Source: 
http: / /www. waterencyclopedia.com/] 

Of the proposed techniques for abating anthropogenic carbon-dioxide, storage in deep geologic 
formations is the only method that has gained widespread acceptance for its feasibility [2|, [3]. 
If geological carbon-dioxide sequestration is implemented on the scale needed to make noticeable 
reductions in atmospheric CO2, a billion metric tons or more must be sequestered annually - a 
250-fold increase over the amount sequestered today. Large sedimentary basins are considered 
best suited to sequester such large volumes of CO2 as they have tremendous pore volume and 
connectivity and they are widely distributed [J, [5] . 

The development of energy systems like geological carbon-dioxide sequestration requires the 
understanding and predictive simulations of complex processes in earth systems. In particular, 
securing a large volume of carbon-dioxide will require a solid scientific foundation defining the cou- 
pled hydrologic-geochemical-geomechanical processes that govern the long term fate of CO2. This 
path is becoming more dependent on modeling coupled geomechanical, thermal, fluid and chemi- 
cal processes and the response of natural and engineered systems. Similarly, these models can be 
utilized to simulate oil and gas reservoirs, advanced recovery from tar sands and shales, under- 
ground storage of hydrogen, natural gas and oil, and evaluating aquifers for new water resources. 
Significant progress in these areas requires new coupled modeling approaches. 

Unfortunately, there are still a number of unanswered scientific questions regarding deep geo- 
logical sequestration that critically affect the efficacy and safety of this method. Although a number 
of studies have been conducted on geochemical interactions within the storage reservoir, and on the 
flow characteristics of the carbon-dioxide plume, relatively little effort has been devoted towards the 

reservoir's structural integrity. The seal created by the cap rock is one of the primary mechanisms 
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that prevents injected carbon-dioxide from escaping back into the atmosphere. Failure of this seal 
can release large quantities of carbon-dioxide, which will have dire consequences. For example, 
in 1986 more than 1700 people were killed by the sudden release of geologic carbon-dioxide from 
lake Nyos in Cameroon (see Figure [T]). Although, in the case of lake Nyos, the carbon-dioxide 
was sequestered naturally, one can learn valuable lessons from this disaster. What happened in 
lake Nyos is similar to opening a shaken bottle of carbonated soda. As long as the cap is on, the 
carbon-dioxide gas stays dissolved under pressure. But when the cap is removed, the bubbles (and 
the soda) rapidly flow out of the bottle. The lake Nyos incident clearly highlights the importance 
of the structural integrity of the cap rock, and an immediate need for a systematic study along the 
lines presented in this paper. 

It is important to note that very few parallel processing commercial and/or research software 
tools exist for simulating complex processes such as coupled multiphase flow with chemical transport 
and geomechanics. Current computational limitations place significant restrictions on realistic 
problems that can be solved. Predictive computational simulation of the coupled-physics associated 
with geosystems is a critical enabling technology for their optimal management. A major stumbling 
block to high-fidelity modeling and simulation of geosystems is the lack of viable, robust, and 
efficient computational technologies for describing multiphase, multicomponent chemically reactive 
fluid mixtures in heterogeneous deformable geologic media. The present paper aims at advancing 
mathematical and numerical modeling, and predictive simulation of flow in deformable porous 
media under high pressures. 

1.1. Main contributions of this paper. Some of the main contributions of this paper are 
as follows: 

(a) We have presented a mathematical model based on the theory of interacting continua which is 
capable of capturing surface subsidence and consolidation of soils. The model is fully coupled, 
and the interaction between the fluid and the porous solid is modeled through drag-like term 
and deformation-dependent porosity. 

(b) We have presented various numerical coupling algorithms, which can be used to obtain the 
coupled response. 

(c) We have presented various representative numerical examples illustrating the predictive capa- 
bilities of the proposed mathematical model, and the numerical performance of the coupling 
algorithms. 

1.2. Organization of the paper. The remainder of this paper is organized as follows. In 
Section[2]we give a brief review of the theory of interacting continua. In Section[3]we present a model 
for the flow of an incompressible fluid in deformable porous solid. In Section U] we shall present 
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Figure 2. Elementary volume of the porous material SV, and the volume occupied by the 
i-th constituent 6V'^'^ (i = 1, • • • , iV). 



coupling algorithms that will be employed to obtained the coupled flow-deformation response. 
Several representative numerical results will be presented in Section [U and conclusions are drawn 
in Section [6l 

2. ELEMENTS OF THE THEORY OF INTERACTING CONTINUA 

We shall employ the mathematical framework offered by the theory of interacting continua 
(TIC), which is sometimes referred to as the theory of mixtures. The basic assumption of the 
theory of interacting continua is that the constituents can be homogenized and assumed to co- 
occupy the domain occupied by the mixture. The theory of interacting continua traces its origins 
to the works of Darcy [B] (also see its English translation by Patricia Bobeck [7]) and Fick [5]. 
Truesdell later gave the theory a firm mathematical footing (see [9], I10|, 111] and several appendices 
in Reference |12j ) . We now provide a brief review of the basic equations of the theory of interacting 
continua. A more detailed treatment can be found in Atkin and Craine |13j . Bowen |14] . and 
Bedford and Drumheller |15j . One can also consult texts by Bowen |16j , Coussy |17|, I18j , de Boer 
|19j . Rajagopal and Tao [20| . and Voyiadjis and Song |21j . 

In the remainder of the paper, it should be noted that the usual summation convention on 
repeated indices will not be adopted. Consider a mixture of N components. One of the components 
can be the porous solid, which will be the case when dealing with a deformable porous solid. A 
typical particle belonging to each constituent in the reference state is denoted X^*^ (i = 1, • • • , N). 
At time t, these particles occupy the position x. The motion and velocity of each constituent 
(z = 1, • • • , N) are defined through 



V 



X 



X«(X«,t) 
dt 



(1) 



(2) 
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We denote the gradient and divergence operators with respect to x through grad[-] and div[-], 
respectively. We now define the volume fractions and porosity by considering a representative 
volume element (RVE) with volume 6V. The volume occupied by the i-th component in the RVE 
is denoted SV^'^h The volume fraction of the i-th component is then defined as follows: 

:= ^ (3) 

By definition, the volume fractions satisfy the following relationship: 

N 



J^^«(^) = l (4) 



i=l 

The porosity of the porous medium is defined as follows: 

cj){x) := 1 - ij^'\x) (5) 

where rj^'^^x) is the volume fraction of the porous solid. Let the mass of the i-th component in this 
RVE be 6m^^^ . The true and bulk mass densities of the i-th component are, respectively, denoted 
by 7^*) and p^'^h That is, 

^ (T) 
The bulk mass density of the mixture is defined as follows: 

N 

p{x,t):=Y,P^^{^.t) (8) 

i=l 

The mixture velocity is defined as follows: 

1 ^ 

^(a.,t):=-^5^p«(^,t).;»(^,t) (9) 

The diffusion velocity v^'^\x,t) of the i-th constituent is defined as follows: 

v^'\x,t) =v^\x,t) -v{x,t) (10) 

Let a{x,t) be any quantity (scalar, vector, tensor) defined at the point x in the mixture at time t. 
We then define 

:,t) +grad[a(a;,t)] • v^''>{x,t) (11a) 
•.,t) + gra,d[a{x^t)] ■ v{x,t) (lib) 





d 


m 




Da 


d 




■~ dt 



The velocity gradient for the i-th constituent L^^^ and its symmetric part Z?^*^ are, respectively, 
defined through 

Z» := grad[^;(*)] (12) 

(13) 

We assume the existence of a partial traction vector and a partial stress tensor T^*^ associated 
with each constituent of the mixture such that 

^(i) = (T»)^n5 (14) 
where ns is the normal to the surface S. We define the total traction and total stress through 

N 

t = ^t« (15) 

i=l 

N 

T = ^T(*) (16) 



i=l 

SO that we have 

t = T^ns (17) 

2.1. Balance laws for constituents and the mixture. The balance of mass for the i-th 
constituent takes the following form: 

^ + divfp^^^-yW] = m(^) (18) 
ot 

where m^*) is the mass supply for the i-th constituent. It is noteworthy that the theory of interacting 
continua can take into account the possible inter-conversion of mass due to chemical reactions 
between the constituents. The balance of mass for the mixture as a whole warrants that 

N 

^m»=0 (19) 

i=l 

The balance of linear momentum of individual constituents can be written as follows: 

^(,) = div[T«]^ + + i» (20) 

where b^^^ is the external body force that acts on the i-th constituent, and i^^^ represents the 
interactive forces acting on the i-th constituent. That is, i^*^ arises due to the forces exerted by the 
other constituents on the i-th constituent by virtue of their being forced to co-occupy the domain 
of the mixture. Constitutive relations need to be specified for these interaction forces and they in 
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general depend on not just the i-th constituent but also on the other constituents. Also, it follows 
from Newton's third law that 



N 

=0 (21) 

1=1 

In the case of a mixture, even in the absence of angular momentum supply, the individual partial 
stresses need not be symmetric. However, the total stress will be symmetric in the absence of 
angular momentum supply, i.e., 

TV /TV \ 

^t(^)= (22) 



1=1 \i=l / 

We shall assume that the individual partial stresses are symmetric, and hence equation (j22p will 
be satisfied automatically. That is, 

rp{t) ^ frp{i)Y Vi =!,•••, TV (23) 



One can similarly write governing equations for the balance of energy and the second law of thermo- 
dynamics. Herein, we shall restrict the studies to isothermal processes, and hence these equations 
are either trivially satisfied or not coupled with the mechanical aspects of the problem. 



3. A MODEL FOR FLUID FLOW IN DEFORMABLE POROUS SOLID 

The above framework offered by the theory of interacting continua is quite general and can 
model a wide variety of interesting phenomena. Herein, however, we shall consider a simplified but 
representative model for fluid and solid, which will be a special case of the above comprehensive 
framework. This simplified model will be used in our study on the stability of coupling algorithms. 
It should be noted that even this simplified model can capture many interesting features in the fiow 
of fluids in deformable solids as discussed later in this paper. 

We do not consider chemical reactions, and hence can set 

=0 Vi = !,••• ,iV (24) 

The constitutive equation for the fluid is taken to be that of a perfect fluid. That is, 

= -p(f)l (25) 

where T^-f) is the stress in the fluid, / is the second-order identity tensor, and p^-^^ is the pressure 
in the fluid. We take into account interactions at the fluid-solid interface by including the following 

7 



drag-like terms: 

iif) ^ ^ifs) f^^if) _ ^(s)^ (26a) 
iis) ^ ^{sf) / (.) _ ^(/)\ (26b) 



By Newton's third law (j2ip . we have The drag term models the friction between 

the fluid and the porous solid, and it does not consider the friction within the fluid. The drag 
coefficient a^'^'*-* can be written as follows: 

a(/^) = ^ (27) 

where is the coefficient of viscosity of the fluid, and k is permeability. In this paper, we shall 
allow the coefficient of viscosity is allowed to depend on the pressure, and is modeled using Barus 
formula [22^ (23]: 



(28) 



where /3*-'^^ has units of Pa^"*^. The values of Hq^"^ and (3^-^^ of various organic liquids can found in 
references [24^ I22j . It is straightforward to show that the balance of linear momentum, given in 
equation (pO|) . becomes 

pif) + grad[^(^)]^;(^) j + a^-^") (v'^^^ - v'^'^^ + grad[p(-^)] = p^^^b^^^ (29) 

Remark 3.1. A simple model that takes into account the friction between the layers of the fluid 
has been proposed by Brinkman j25j . The model is referred to as Brinkman equation or sometimes 
Darcy-Stokes equation. Under this model, the partial stress in the fluid takes the following form: 

rpif) ^ _p[f)j + 2/i(-^)D(-^) (30) 



In such a case, the above expression for balance of linear momentum (j29p should be altered as 
follows: 

p^f^ + grad[i;(-^V^-^^^ + (^^^^ - v^'^) + gradb^^)] - divp^^-^^D^-^)] = p^^^h^^^ (31) 

The fluid is assumed to be incompressible, and hence the balance of mass for the fluid can be 
written as follows: 



(32) 



where (p is the porosity of the porous media. It should be noted that the formulation can be 
easily extended to variable density flows. We shall model the subsurface rock as a linearized elastic 



material. That is, 



max ||grad[u('-']||^ < 1 (33) 



where || • ||^, denotes the trace norm [26], which is defined as follows: 

Va^a 



tr 



(34) 



where A is a second-order tensor, and tr[-] denotes the trace. The (linearized) strain in the porous 
solid is given by 

e(^) := ^ (grad[u('*)] + grad[-u(^)]^) (35) 
Note that the velocity of the solid is given by 

The governing equation for the balance of linear momentum for the porous solid takes the following 
form: 



9*2 



^W^^ _ «(/^) (^(/) _ ^(^)) _ diviT^'')] = p(^)b('^) (37) 



where u^^^ is the displacement of the solid. Note that the interaction effect has opposite signs in 
the two balance equations (p9|) and ([37|) in virtue of Newton's third law. The constitutive equation 
for the porous solid can be written as follows: 

T(^) = A(^Hr[e(^)]/ + 2/iWe(") (38) 

where T^*-* is the stress in the porous solid, and A^*^ and [i^^^ are the Lame parameters of the porous 
solid. A few remarks about the above mathematical model are in order. 

Remark 3.2. The partial stresses for the fluid (given by equation ()25p ) and for the solid (given 
by equation (j38p ) are symmetric, which is in accordance with the assumption discussed in the 
previous section. 

Remark 3.3. Equations (I29p . (I32p and (I37p are coupled and are in terms of three unknowns: 
p(f) ^ v^f\ and u^'^\ In some of the specialized models discussed below, the number of unknowns can 
be reduced by substituting one or more of the equations in the above system into the other equations. 



3.1. Quasi-static models. One popular approach to incorporate the influence of the fluid 
motion on the deformation of the solid and vice versa is by modeling the porosity of the porous 
media to be dependent on the deformation gradient in the solid. Typically, in this scenario, a 
number of approximations are made which are necessary to model the balance of momentum for 
the fluid using Darcy's law. If we ignore inertia in both the fluid and the solid, we can drop the first 
term in equations (j29p and (j37p . Similarly, if we assume that the solid velocity is much smaller than 
the fluid velocity, the interaction term (second term in equations (j29p and (|37p) may be written as 

Darcy's equation can be substituted into the interaction term for the fluid velocity. Rather 
than replace the interaction, —a^f^'^v^f\ with the gradient of the fluid pressure, grad[p^'^^], and the 
fluid body force, —p^f^b^-f \ we shall incorporate the influence of the fluid pressure in the solid stress 
tensor as 

where Tp^ = —p^^^I. Partially incorporating the interaction term due to the fluid pressure in the 
solid stress as above is possible due to the following relationship, grad[p'^''^^] = — div[Tp''^]. The fluid 
body force, p^f^h^f \ must of course be taken into account. 

3.2. Steady-state response. The governing equations for steady-state response of flow of an 
incompressible fluid in a deformable porous solid can be written as follows: 

Governing equations for fluid 

a{/«)(p{/))^(/) + gradb^^)] = (40a) 

div[0'y(-^)] = (40b) 
Governing equations for solid 

div[T(^)] + a^f'\p^f^)v'^f^ + = (40c) 

T(") = A("Hr[e(")] J + 2^(")e('^) (40d) 

^(s) ._ 1 /'grad[it(^)] + grad[it(^)]^) (40e) 



2 

Deformation-dependent porosity 



(40f) 



1 + (1 - (l)oMe(^^] 

where (I)q is the porosity when the porous solid is unstrained (i.e., tr[e(*)] = 0). The unknown field 
variables are: velocity of the fiuid v^'f^x), pressure in the fluid p^f\x), displacement of the porous 
solid u^^\x), and porosity <^(a;). The above model (given by equations (I40ap ~ (|40fp ) is nonlinear 
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and fully coupled in the sense that neither the fluid model nor the solid model can be solved 
independently. 

Remark 3.4. In the literature, one can find many other models for deformation- dependent 
porosity. For example, the following expression has been proposed: 

exp[tr[e^''-'JJ 

Another model, which is used for large deformations of the porous solid, takes the following form: 

^ = cJl-^^^] (42) 

where (po is the initial porosity, F^^^ is the solid deformation gradient 

pis) ^ J ^ grad[it(^)] (43) 

and C(j) is a coefficient that may very in time, and can be used to smooth jumps in the porosity in 
time. In the case where the mass balance equation ( equation ()40bp ) is subcycled, large changes in 
porosity can occur during synchronization updates. Using a reference pressure, p^d; that is taken as 
the pressure at the last synchronization update, the following expression can be used to interpolate 
the porosity between updates 

C^ = l + Cr{p-p,,i) (44) 

where Cr is the rock compressibility constant. 

4. NUMERICAL COUPLING ALGORITHMS 

The coupled equations presented above can be solved in a number of ways, which will depend on 

how the overall problem is decomposed. The first method is the fully coupled method, which solves 

governing equations for fluid, solid and porosity update simultaneously as one monolithic system. 

The variables in this case are solid displacement, fluid velocity, fluid pressure, and porosity. Under 

this method, it is typical to use the same time step for all subsystems. The second method is the 

lockstep method (which is also referred to as Gauss-Seidel), which solves each subsystem individually 

until convergence in serial fashion. The solution is then transferred to the other subsystems which 

then iterate until convergence. Under the lockstep method, the same time step is also used for all 

subsystems. In the third method, referred to as the subcycle method, subcycling occurs in either 

the mass balance or solid deformation equation. For the problems presented below, subcycling was 

only used for the mass balance equation since the physics of the flow evolves much faster than the 

solid deformation. The last method considered in this paper is the Jacobi method. Under this 

method, both equations are decoupled and converged individually, but fields are transferred after 

each iteration. Figure [3] pictorially describes these coupling algorithms. 
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5. REPRESENTATIVE NUMERICAL RESULTS 

In this section, using the mathematical model described in Section [Sj we shall study the perfor- 
mance and relative efficiency of the coupling algorithms that are outlined in the previous section. 
Several canonical problems are solved. 

5.1. Verification by a manufactured solution. In this subsection, we shall use the method 
of manufactured solutions to verify the implementation of the coupling algorithms. The method of 
manufactured solutions assumes an exact solution and then derives the corresponding source terms 
that satisfy the governing equations. The (manufactured) exact solution for the coupled problem 
(|40ap - (|40fp takes the following form: 

= sin(7rx), €^'\x) = u^'K cos{ttx), T^'\x) = + 2^(^)) 4'^7r cos(7rx) (45a) 

v^f\x) = v^J'^ (l + (1 - </'o)4'^vrcos(^x)) , p^^^x) = -a^f'^v^/'^ [x + {1 - (/.o)4'^ sin(^x)) + p^f^ 

(45b) 

6(x) = ^ = (45c) 

1 + (1 - <^o)tr[e(-)] 1 + (1 _ 0o)4'^7r cos(^x) 

The prescribed body forces for the fluid and the porous solid are as follows: 

(s) 2 

b^f\x) = 0, b^'\x) = ^^^j^ (a^**) + 2;uW) sin(7rx) (46) 
The constants used in this numerical experiment are listed in Table [TJ To restrict the analysis of 
Table 1. Verification by a manufactured solution: Values used in the numerical simulation. 



Parameter Value Parameter Value 



pif) 


1.0 


pis) 


1.0 


^0 


1.0 


„(/) 


1.0 


A(-) 


1.0 




0.5 


«(/^) 


1.0 


4>o 


0.1 



lif^ 0.01 



this problem to one-dimension, the boundary conditions are prescribed such that displacements in 
the y-direction are fixed and only displacements in the x-direction are considered. The problem 
domain consists of a unit square with a single element in the y-direction and 200 elements in the 
x-direction. On the left side of the domain the solid displacement is fixed. On the right side of the 
domain, the manufactured solution is applied for the solid displacements as a boundary condition. 

The above outlined problem was solved by each of the coupling algorithms described in Figure[3j 
A comparison between the exact solution and the obtained numerical solution for the fluid pressure, 
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solid displacements, and the porosity are shown in Figure [H Figures [5] and [6] provide a measure 
of the efficiency of the loosely coupled algorithms (lockstep and Jacobi). The relative error of the 
solution is computed as the L2 norm of the exact solution minus the computed solution. Notice 
that both the lockstep and Jacobi methods perform similarly in terms of the accuracy obtained for 
a given convergence tolerance, but the Jacobi method requires more iterations for a given level of 
accuracy. 

5.2. Terzaghi consolidation problem. To evaluate the performance of the quasi-static 
model (equations ()40ap to ()40fp ) for one-dimensional transient wave propagation in fluid satu- 
rated incompressible porous media we investigate the response of an infinite half-space to time 
dependent loading. This example problem, which has been studied by a number of researchers 
including [27(, I28| . is represented computationally by the domain and parameters shown in Figure 
[71 although the dimensions of the domain are insignificant regarding the resulting response. The 
top surface of the domain is subjected to a forcing function, F, given as: 

F = 100((1 - cos(75t)) (47) 

The distributed force is applied in a weak fashion over the top surface. This surface is also treated 
as perfectly drained meaning that the fluid is free to drain from the top surface and the pressure is 
equal to the ambient pressure. Along the sides and bottom of the domain the fluid is not allowed to 
drain. The displacement boundary conditions consist of fixing the x-displacement along the right 
and left boundary and the y-displacement on the bottom surface. This essentially allows for only 
one-dimensional consolidation in the y-direction. 

The analytic solution for the y-displacement of the top surface of the domain is given as: 

w = vs(A.:."2,..)) i* - ^>'"^'» i^"^^) - <^«' 

where a = (n(^)2p^^^ +n(/)2pW)/c, h = n^f'i'^kjc, c = (A^ + 2/iW)n(/)2, = I x 10^ Iq{z) is 
the modified Bessel function of zeroth order, and U{t) is a unit Heaviside function. The computed 
solution using the proposed quasi-static model is shown in Figure[51 Notice that the model performs 
well even though the loading is nonlinear and the dynamic response of the solid is treated in a quasi- 
static fashion. 

5.3. Surface subsidence problem. In this verification problem we consider the subsidence 
of a reservoir resulting from extraction of fluid via a centrally located well that has been studied 
extensively in [29j . Figure O shows the domain and boundary conditions used for this problem. In 
order to create confinement conditions similar to the subsurface environment an initial stress state 
was specified for the solid stress. The coupling between the solid stress and the fiuid pore pressure 
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is defined by equation (j39p and the porosity model corresponds to equation ()40f|) . The quasi-static 
model given in the previous examples was also used for the governing equations. Figure [TU] shows a 
comparison between the results published in j29j and the computed results from the present study. 
Results are shown for both coupled and uncoupled porosity. Close correlation is achieved for the 
subsidence of the center of the top surface the domain for the coupled results. If the effect of the 
solid deformation is not taken into account in the porosity, the subsidence of the top surface is 
under-predicted by roughly 5%. Figure [11] shows contours of porosity for the top surface of the 
domain. Notice that near the ejection well, the porosity is decreased by almost 30% from the initial 
value. 

5.4. Two-phase immiscible water flood. In the last example we present results for a water 
flood of the five spot well domain. The problem domain consists of two-dimensional square with 
injection wells at one of the sets of opposing corners and ejection wells at the other. Water is 
injected into the domain pushing oil out the ejection wells. The solution to this problem has been 
published for a number of numerical approaches including the results found in |30j . 

5.4.1. Changes in porosity vs. change in permeability due to solid stress. In the present study 
we consider the effects of both changes in porosity due to stresses in the solid skeleton and resulting 
changes in permeability due to damage. Conceptually, it is well established that solid deformation 
dependent porosity models behave like compressibility models. The otherwise incompressible fluid 
is able to be squeezed into pores that are expanding and vice versa. This is also apparent from the 
solid stress dependent contribution to the fluid equations showing up in the mass balance equation. 
In addition to the effect of pores enlarging and collapsing, one may wish to incorporate the increased 
permeability or interconnectedness of the void spaces due to solid damage. For example, in regions 
of high stress in the reservoir rock, cracks will form that allow fluid to flow with much less resistance. 
As a result, the flow magnitude and path of the fluid may change depending on the state of stress in 
the geologic formation. To explore differences between porosity coupled and permeability coupled 
models we investigate the two-phase water flood problem using both models separately. For the 
porosity coupled results we employ the model given in equation ()40f p . for the permeability coupled 
results we consider a damage modified permeability model given as: 



where is a scaling parameter. The permeability coupled model is based on the magnitude of the 



situ conditions, ||Tq ||. This ad-hoc model will serve as the basis for comparison with the coupled 
porosity model. In a forthcoming work, we investigate damage modified permeability models more 
thoroughly. 




(49) 
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Since we are dealing with two-phase immiscible flow, we modify the governing equations of the 
quasi-static model as follows: 



(p(/)<A(l-S„))+div[p(/)^(/)] 


= in 


(50a) 


|(pi/)<^5„) + div[p(/)t,(/)] 


= m^^ in Q 


(50b) 


div[T(^)] 


= inn 


(50c) 




= pP{x) on TP 


(50d) 




= uP{x) onT" 


(50e) 




= r onr* 


(50f) 



where S is the saturation and the subscripts "w" and "n" denote the wetting and non-wetting 
phases respectively. Note that we have incorporated the relationship Sn + S^] = 1.0 in the first 
equation. The Darcy velocities for this problem are given as: 

v[i^ = -«(/)(grad[p(/)] + p(/)6(/)) (51) 
^(/) = -a(-«(grad[p(/)] + ^^if^) (52) 

Here we have assumed no capillary pressure such that pn = Pw With this assumption, the equations 
above may be solved for Sn and pn (or pw). Note that we are using the same damage modified 
permeability a^-^-* for both phases. To limit the spurious oscillations engendered by the advection 
terms, we use the CVFEM upwinding method presented by Forsyth in references [311, 132| . 

The problem domain and boundary conditions are shown in Figure [T2l The elastic properties for 
the reservoir (used only in the coupled results) were generated by the code Sierra/Encore |33L I34j 
using a Karhunen-Loeve realization. Contours of A^*-* are shown in Figure [THl Notice large regions 
of variation of the elastic properties throughout the domain. Also shown in Figure [13] are contours 
of the porosity and permeability at time t = 2 x 10^ s for both porosity coupled and permeability 
coupled results. Figure [T3] shows the tensor magnitude contours of the solid stress due to the 
pore pressure loading. When compared with Figure [13] it reveals that the porosity model exhibits 
its largest porosity in a concentrated region near the region of low solid stress magnitude. The 
permeability model on the other hand shows high permeability in regions of high stress as the 
model intends. This results in a drastically different flow behavior between the models as shown 
in Figure [15] The top-left side of Figure [15] shows the non- wetting phase saturation for uncoupled 
two phase model. The top-right part of the figure shows the non-wetting phase saturation for the 
porosity coupled model. Notice that the differences with the uncoupled model are negligible. In 
contrast, the damage modified permeability model results clear show race-tracking in the regions 

15 



of high stress, which is shown the bottom side of Figure [151 This example highhghts the need for 
further model development to capture damage motivated permeability effects. 

6. CONCLUDING REMARKS 

We have considered a comprehensive mathematical model for flow of an incompressible fluid in 
deformable porous solid. The model is based on the theory of interacting continua with deformation- 
dependent porosity and damage modified permeability. The model also takes into account the 
dependence of viscosity the fiuid on the pressure of the fiuid. The model is fully coupled in the 
sense that the flow problem depends on the response of the porous solid, and the deformation of the 
porous solid depends on the velocity of the fluid. Several numerical coupling algorithms to obtain 
coupled flow-deformation response are also presented. Representative numerical results show the 
importance of considering the deformation of the solid on the flow and vice- versa. A possible future 
work is develop models and numerical coupling algorithms for flow in an elasto-plastic fracturable 
porous solid. 

ACKNOWLEDGMENTS 

This work was supported in part by Sandia National Laboratories through the Laboratory 
Directed Research and Development program (Contract No. Cll-00239). Sandia is a multipro- 
gram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United 
States Department of Energy's National Nuclear Security Administration under Contract DE- 
AC04-94AL85000. The second author (K. B. Nakshatrala) also acknowledges the support of the 
National Science Foundation under Grant no. CMMI 1068181. The opinions expressed in this 
paper are those of the authors and do not necessarily reflect that of the sponsors. Dr. Martinez was 
supported in part by the Center for Frontiers of Subsurface Energy Security, an Energy Frontier 
Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy 
Sciences under Award Number DESC0001114. 

References 

[1] A. Leach, C. F. Mason, and K. Veld. Co-optimization of enhanced oil recovery and carbon sequestration. Tech- 
nical Report NBER Working Paper No. 15035, National Bureau of Economic Research, 2009. 
[2] A. M. Macfarlance. Energy: The issue of the 21''' century. Elements, 3:165-170, 2007. 
[3] D. P. Schrag. Confronting the climate-energy challenge. Elements, 3:171-178, 2007. 

[4] S. Bachu. Screening and ranking of sedimentary basins for sequestration of CO2 in geological media in response 

to climate change. Environmental Geology, 44:277-289, 2003. 
[5] S. M. Benson and D. R. Cole. CO2 sequestration in deep sedimentary formations. Elements, 4:325-331, 2008. 
[6] H. Darcy. Les Fontaines Publiques de la Ville de Dijon. Victor Dalmont, Paris, 1856. 

[7] P. Bobeck. The Public Fountains of the City of Dijon. Kendall Hunt Publishing, Dubuque, Iowa, USA, 2004. 

16 



[8] A. Fick. Uber diffusion. Poggendorff's Annalen der Physik und Chemie, 94:59-86, 1855. 
[9] C. Truesdell. Sulla basi della thermomechanica. Rendiconti dei Lincei, 22:33-38, 1957. 
[10] C. Truesdell. Sulla basi della thermomechanica. Rendiconti dei Lincei, 22:158-166, 1957. 

[11] C. Truesdell and R. A. Toupin. The Classical Field Theories in Handbuch der Physik, volume III/l, pages 

226-793. Springer- Verlag, New York, 1960. 
[12] C. Truesdell. Rational Thermodynamics. Springer, New York, USA, 1984. 

[13] R. J. Atkin and R. E. Craine. Continuum theories of mixtures: Basic theory and historical development. The 

Quaterly Journal of Mechanics and Applied Mathematics, 29:209-244, 1976. 
[14] R. M. Bowen. Theory of mixtures. In A. C. Eringen, editor. Continuum Physics, volume III. Academic Press, 

New York, 1976. 

[15] A. Bedford and D. S. Drumheller. Theories of immiscible and structured mixtures. International Journal of 

Engineering Science, 21:863-960, 1983. 
[16] R. Bowen. Porous Elasticity: Lectures on the Elasticity of Porous Materials as an Application of the Theory of 

Mixtures. Available online: http://repository.tamu.edu/handle/1969. 1/2500, Texas A&M University, 2010. 
[17] O. Coussy. Poromechanics. John Wiley & Sons, Inc., New York, USA, 2004. 

[18] O. Coussy. Mechanics and Physics of Porous Solids. John Wiley & Sons, Inc., New York, USA, 2010. 

[19] R. de Boer. Theory of Porous Media: Highlights in the Historical Development and Current State. Springer- 

Verlag, New York, USA, 2000. 
[20] K. R. Rajagopal and L. Tao. Mechanics of Mixtures. World Scientific Publishing, River Edge, New Jersey, USA, 

1995. 

[21] G. Z. Voyiadjis and C. R. Song. The Coupled Theory of Mixtures in Geomechanics with Applications. Springer- 

Verlag, New York, USA, 2006. 
[22] C. Barus. Isotherms, isopiestics and isometrics relative to viscosity. American Journal of Science, 45:87-96, 1893. 
[23] K. B. Nakshatrala and K. R. Rajagopal. A numerical study of fluids with pressure dependent viscosity flowing 

through a rigid porous medium. International Journal for Numerical Methods in Fluids, 67:342-368, 2011. 
[24] P. W. Bridgman. The Physics of High Pressure. MacMillan Company, New York, USA, 1931. 
[25] H. C. Brinkman. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. 

Applied Scientific Research, Al:27-34, 1947. 
[26] P. D. Lax. Functional Analysis. John Wiley & Sons, Inc., New York, USA, 2002. 

[27] R. de Boer, W. Ehlers, and Z. Liu. One-dimensional transient wave propagation in fluid-saturated incompressible 

porous media. Applied Mechanics, 63:59-72, 1993. 
[28] S. Diebels and W. Ehlers. Dynamic analysis of a fully saturated porous medium accounting for geometrical and 

material non-linearities. International Journal for Numerical Methods in Engineering, 39:81-97, 1996. 
[29] R. H. Dean, X. Gai, C. M. Stone, and S. E. Minkoff. A comparison of techniques for coupling porous flow and 

geomechanics. Society of Petroleum Engineers, SPE 79709:132-140, 2003. 
[30] R. Helmig. Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of 

Hydrosystems. Springer- Verlag, Berlin, second edition, 1997. 
[31] P. Forsyth. A control-volume, finite-element method for local mesh refinement in thermal reservoir simulation. 

SPE Reservoir Engineering, 5:561-566, 1990. 
[32] P. Forsyth. A control volume finite element approach to NAPL groundwater contamination. SIAM Journal on 

Scientific and Statistical Computing, 12:1029-1057, 1991. 



17 



[33] P. Notz, S. R. Subia, M. Hopkins, H. K. MofTat, and D. R. Noble. Aria 1.5 User Manual. Technical report, Sandia 

Report SAND2007-2734, Sandia National Laboratories, 2007. 
[34] Fuego. SIERRA/Fuego 2.7 User's Manual. Sandia National Laboratories, SAND 2006-6084P, 2008. 



18 



Dr. Daniel Z. Turner, Department of Civil Engineering, University of Stellenbosch, Stellen- 
BOSCH, South Africa. TEL: +27-21-808-4434 
E-mail address: dzturner@sun.ac.za 

Dr. Kalyana Babu Nakshatrala, Department of Civil and Environmental Engineering, University 
OF Houston, Houston, Texas - 77204. TEL: + 1-713- 743-4418 
E-mail address: knakshatrala@uh.edu 

Dr. Mario J. Martinez, Thermal and Fluid Processes Department, Sandia National Laboratories, 
Albuquerque, New Mexico - 87185-0836. TEL: + l-505-844-8729 
E-mail address: mjmarti@sandia.gov 



19 





-| Time step J 1 






Solve fluid 
and solid 


















Figure 3. Numerical coupling algorithms: fully coupled (top left), lockstep (top right), subcycling (bottom left), and Jacobi 
(bottom right). 
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Figure 4. Verification by a manufactured solution: numerical solution vs. x compared 
with the exact solution (top left) solid displacement, Ux (top right) pressure, p (bottom 
left) fluid velocity, Vx (bottom right) porosity, (j). The lockstep method was computed using 
a tolerance of e; = 1 x 10~^. The subcycle method was computed using one subcycle 
(equivalent to lockstep) and a tolerance of Es = 1 x 10~^. The Jacobi method was computed 
using a tolerance of = 1 x 10~^. 
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Figure 5. Verification by a manufactured soiution: minimum error vs. global convergence 
tolerance, e; or ej . The error is measured as the i2-norm of the computed solution minus 
the exact solution and the minimum is taken over the total solution time. 
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Figure 6. Verification by a manufactured solution: number of iterations vs. global conver- 
gence tolerance, e; or Cj. The number of iterations is measured as the number of fluid and 
solid solves per solution cycle before the global tolerance is met. 
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Figure 7. Terzaghi dynamic consolidation: problem geometry and boundary conditions 




Figure 8. Terzaghi dynamic consolidation: computed solution compared with the analytic 
solution given in 
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Figure 9. Surface subsidence: Problem geometry, parameters, and boundary conditions. 
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Figure 10. Surface subsidence: Comparison of the subsidence at the center of the top of 
the domain with pubhshed results in [29| for both coupled and uncoupled porosity. 
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Figure 11. Surface subsidence: Contours of porosity for the top surface of the domain. 



,(/) 



M^f' = 0.0 
wlf' = 0.0 




m. 



(/) 



,(/) 





(mean) 


= 5.9402 X 10^ 




(mean) 


= 5.767 X 10** 




pU) 


= 1000 






= 1 X 10(-^) 




00 


= 0.2 






= 1.0 




P 


= 2.0 






= 3.2 X 10"^ 






= -3.2 X 10"^ 



Figure 12. Two phase immiscible water flood: Problem geometry, parameters, and boimd- 
ary conditions. 
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Figure 13. Two phase immiscible water flood: 
and 4> (bottom). 
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Figure 15. Two phase immiscible water flood: Contom-s of Sn- The top-left figure shows 
the non- wetting phase saturation for uncoupled two phase model. The top-right figure shows 
the non-wetting phase saturation for the porosity coupled model. The bottom figure shows 
the non-wetting phase saturation obtained using the damage modified permeability model. 
There are no noticeable differences between the results between the uncoupled model and 
the porosity coupled model. In contrast, the damage modified permeability model results 
clear show race-tracking in the regions of high stress. 
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